function [rSVtTEItR, Ek, vSVtTEItR, EkRate] = svposvel_gps(tR, tau, ephemPar, cst, nEccAnomIter) %#codegen
% SVPOSVEL_GPS 根据星历计算导航卫星在信号接收的位置和信号发射时刻的速度
%
% Description
% 根据接收机给出的导航电文中的星历计算信号发射时刻的卫星位置、速度在信号接收时刻的ECEF坐标系中的分量
%
% Tips
% # 该函数计算得到的卫星轨道三维位置误差的均方差一般在3~5米以内
%
% Input Arguments
% # tR: 标量，信号接收时刻，单位：s
% # tau: 标量，信号传播时延，单位：s
% # ephemPar: 结构体，包含以下域
%     tOE: 标量，星历数据参考时刻，值域：0~604784，单位：s
%     ASqrt: 标量，长半轴的平方根，单位：m^0.5
%     e: 标量，偏心率，值域：0~0.03，单位：无量纲
%     omega: 标量，近地点幅角，单位：rad
%     Deltan: 标量，卫星平均运动速率与计算值之差，单位：rad
%     M0: 标量，参考时刻的平近点角，单位：rad
%     Omega0: 标量，周内时为0时的升交点赤经，单位：rad
%     OmegaRate: 标量，升交点赤经的变化率，单位：rad/s
%     i0: 标量，参考时刻的轨道倾角，单位：rad
%     IDOT: 标量，轨道倾角的变化率，单位：rad/s
%     Cuc: 标量，纬度幅角的余弦调和改正项的振幅，单位：rad
%     Cus: 标量，纬度幅角的正弦调和改正项的振幅，单位：rad
%     Crc: 标量，轨道半径的余弦调和改正项的振幅，单位：m
%     Crs: 标量，轨道半径的正弦调和改正项的振幅，单位：m
%     Cic: 标量，轨道倾角的余弦调和改正项的振幅，单位：rad
%     Cis: 标量，轨道倾角的正弦调和改正项的振幅，单位：rad
% # cst: 结构体，包含以下域
%     mu: 地球（包含大气层）引力常数，单位：m^3/s^2
%     OmegaeRate: 地球自转角速率，单位：rad/s
% # nEccAnomIter: 标量，计算偏近点角Ek时的迭代次数，默认为10
%
% Output Arguments
% # rSVtTEItR: 3*1矢量，信号发射时刻的卫星位置在信号接收时刻的ECEF坐标系中的坐标，单位：m
% # vSVtTEItR: 3*1矢量，信号发射时刻的卫星速度在信号接收时刻的ECEF坐标系中的分量，单位：m/s
% # Ek：标量，信号发射时刻的偏近点角，单位：rad
% # EkRate：标量，信号发射时刻卫星偏近点角对时间的导数，单位：rad/s
%
% Assumptions and Limitations
% # 本函数适用于GPS卫星及北斗MEO/IGSO卫星
% # 星历参数tOE为星历的参考时刻，计算卫星轨道的时刻应该在tOE前后两小时以内，超过这个范围结果无效
%
% References
% #《应用导航算法工程基础》“GPS卫星轨道计算方法”

if nargin < 5
    nEccAnomIter = 10; % 设置迭代次数
end

%% 计算规化时间tk
tk = tR - tau - ephemPar.tOE;
[tk] = accountforweekcrossover(tk);
if abs(tk) >= 7200
    warning('svposvel_gps:abnormaltimeofephem', 'Absolute value of tk is equal or greater than 7200s.');
end
if tk >= 0
    warning('svposvel_gps:abnormaltimeofephem', 'tk is equal or greater than 0s.');
end

%% 计算卫星的位置
% 计算卫星的平均角速率n
A = ephemPar.ASqrt^2;
n0 = sqrt(cst.mu / A^3);
n = n0 + ephemPar.Deltan;

% 计算信号发射时刻的平近点角Mk
Mk = ephemPar.M0 + n*tk;
if Mk < 0
    Mk = Mk + 2*pi;
end

% 计算信号发射时刻的偏近点角Ek（迭代法解超越方程）
Ek = Mk; % 迭代初始值选择为平近点角
for j = 1 : nEccAnomIter
    Ek = Mk + ephemPar.e*sin(Ek);
end

% 计算信号发射时刻的真近点角vk
cosvk = (cos(Ek)-ephemPar.e) / (1-ephemPar.e*cos(Ek));
sinvk = sqrt(1 - ephemPar.e^2)*sin(Ek) / (1-ephemPar.e*cos(Ek));
vk = atan2(sinvk, cosvk);
Ek = acos((ephemPar.e+cosvk) / (1+ephemPar.e*cosvk));

% 计算信号发射时刻的纬度幅角Φk
Phik = vk + ephemPar.omega;

% 计算信号发射时刻的摄动校正项
deltauk = ephemPar.Cus*sin(2*Phik) + ephemPar.Cuc*cos(2*Phik);
deltark = ephemPar.Crs*sin(2*Phik) + ephemPar.Crc*cos(2*Phik);
deltaik = ephemPar.Cis*sin(2*Phik) + ephemPar.Cic*cos(2*Phik);

% 计算摄动校正后的纬度幅角uk、轨道半径rk、轨道倾角ik
uk = Phik + deltauk;
rk = A*(1-ephemPar.e*cos(Ek)) + deltark;
ik = ephemPar.i0 + ephemPar.IDOT*tk + deltaik;

% 计算信号发射时刻卫星在轨道平面的位置(xkPrime, ykPrime)
xkPrime = rk*cos(uk);
ykPrime = rk*sin(uk);

% 计算信号发射时刻的升交点赤经Ωk
Omegak = ephemPar.Omega0 + (ephemPar.OmegaRate-cst.OmegaeRate)*tk - cst.OmegaeRate*ephemPar.tOE;

% 计算信号发射时刻卫星在ECEF坐标系中的坐标(xk, yk, zk)
xk = xkPrime*cos(Omegak) - ykPrime*cos(ik)*sin(Omegak);
yk = xkPrime*sin(Omegak) + ykPrime*cos(ik)*cos(Omegak);
zk = ykPrime*sin(ik);

% 计算信号发射时刻的卫星位置在信号接收时刻的ECEF坐标系中的坐标
theta = cst.OmegaeRate * tau;
CEItTEItR = [cos(theta) sin(theta) 0; -sin(theta) cos(theta) 0; 0 0 1];
rSVtTEItR = CEItTEItR * [xk yk zk]';

if nargout > 2
    %% 计算卫星的速度
    % 计算信号发射时刻卫星平近点角对时间的导数
    MkRate = n;
    
    % 计算信号发射时刻卫星偏近点角对时间的导数
    EkRate = MkRate / (1-ephemPar.e*cos(Ek));
    
    % 计算信号发射时刻卫星真近点角对时间的导数
    vkRate = sqrt(1-ephemPar.e^2)*EkRate / (1-ephemPar.e*cos(Ek));
    
    % 计算信号发射时刻卫星升交点角距对时间的导数
    PhikRate = vkRate;
    
    % 计算信号发射时刻卫星的摄动校正项对时间的导数
    deltaukRate = 2 * PhikRate * (ephemPar.Cus*cos(2*Phik) - ephemPar.Cuc*sin(2*Phik));
    deltarkRate = 2 * PhikRate * (ephemPar.Crs*cos(2*Phik) - ephemPar.Crc*sin(2*Phik));
    deltaikRate = 2 * PhikRate * (ephemPar.Cis*cos(2*Phik) - ephemPar.Cic*sin(2*Phik));
    
    % 计算改正后的纬度幅角、轨道半径和轨道倾角对时间的导数
    ukRate = PhikRate + deltaukRate;
    rkRate = A*ephemPar.e*sin(Ek)*EkRate + deltarkRate;
    ikRate = ephemPar.IDOT + deltaikRate;
    
    % 计算信号发射时刻卫星在轨道平面的位置对时间的导数
    xkPrimeRate = rkRate*cos(uk) - rk*ukRate*sin(uk);
    ykPrimeRate = rkRate*sin(uk) + rk*ukRate*cos(uk);
    
    % 计算信号发射时刻卫星的升交点赤经对时间的导数
    OmegakRate = ephemPar.OmegaRate - cst.OmegaeRate;
    
    % 计算信号发射时刻卫星的速度(vxk,vyk,vzk)
    xkRate = xkPrimeRate*cos(Omegak) - yk*OmegakRate - (ykPrimeRate*cos(ik)-zk*ikRate)*sin(Omegak);
    ykRate = xkPrimeRate*sin(Omegak) + xk*OmegakRate + (ykPrimeRate*cos(ik) - zk*ikRate)*cos(Omegak);
    zkRate = ykPrimeRate*sin(ik) + ykPrime*cos(ik)*ikRate;
    
    % 计算信号发射时刻的卫星速度在信号接收时刻的ECEF坐标系中的分量
    vSVtTEItR = CEItTEItR * [xkRate ykRate zkRate]';
end
end